# defining an operator for "not in"

'%!in%' <- function(x,y)!('%in%'(x,y))

# Function to add a common legend for combined ggplots

g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}

# function to install packages (if necessary) and load all of them

usePackage<-function(p){
  # load a package if installed, else load after installation.
  # Args:
  #   p: package name in quotes
  for(i in unique(p)){
    if (!is.element(i, installed.packages()[,1])){
      print(paste('Package:',i,'Not found, Installing Now...'))
      install.packages(i, dep = TRUE)}
    print(paste('Loading Package :',i))
    require(i, character.only = TRUE) 
  }
}

# modified function from idealstan package to estimate ideal points

estimate_ideal_points <- function (idealdata = NULL, model_type = 2, inflate_zero = FALSE, 
                                   vary_ideal_pts = "none", keep_param = NULL, grainsize = 1, 
                                   mpi_export = NULL, use_subset = FALSE, sample_it = FALSE, 
                                   subset_group = NULL, subset_person = NULL, sample_size = 20, 
                                   nchains = 4, niters = 1000, use_vb = FALSE, ignore_db = NULL, 
                                   restrict_ind_high = NULL, fix_high = 1, fix_low = (-1), restrict_ind_low = NULL, 
                                   fixtype = "prefix", const_type = "persons", id_refresh = 0, 
                                   prior_fit = NULL, warmup = 1000, ncores = 4, use_groups = FALSE, 
                                   discrim_reg_sd = 2, discrim_miss_sd = 2, person_sd = 3, time_fix_sd = 0.1, 
                                   time_var = 10, ar1_up = 1, ar1_down = 0, boundary_prior = NULL, 
                                   time_center_cutoff = 50, restrict_var = FALSE, sample_stationary = FALSE, 
                                   ar_sd = 1, diff_reg_sd = 1, diff_miss_sd = 1, restrict_sd_high = 0.01, 
                                   restrict_sd_low = 0.01, tol_rel_obj = 0.001, gp_sd_par = 0.025, 
                                   gp_num_diff = 3, gp_m_sd_par = 0.3, gp_min_length = 0, cmdstan_path_user = NULL, 
                                   gpu = FALSE, map_over_id = "persons", save_files = NULL, 
                                   pos_discrim = FALSE, het_var = TRUE, compile_optim = TRUE, 
                                   ...) 
{
  if (is.null(cmdstanr::cmdstan_version())) {
    print("You need to install cmdstan with cmdstanr to compile models. Use the function install_cmdstan() in the cmdstanr package.")
  }
  if (use_subset == TRUE || sample_it == TRUE) {
    idealdata <- subset_ideal(idealdata, use_subset = use_subset, 
                              sample_it = sample_it, subset_group = subset_group, 
                              subset_person = subset_person, sample_size = sample_size)
  }
  if (!is.null(cmdstan_path_user)) 
    set_cmdstan_path(cmdstan_path_user)
  if (!file.exists(system.file("stan_files", "irt_standard_map", 
                               package = "idealstan"))) {
    print("Compiling model. Will take some time as this is the first time the package has been used.")
    print("Have you thought about donating to relief for victims of Yemen's famine?")
    print("Check out https://www.unicef.org/emergencies/yemen-crisis for more info.")
  }
  stan_code_map <- system.file("stan_files", "irt_standard_map.stan", 
                               package = "idealstan")
  stan_code_gpu <- system.file("stan_files", "irt_standard_gpu.stan", 
                               package = "idealstan")
  if (compile_optim) {
    idealdata@stanmodel_map <- stan_code_map %>% cmdstan_model(include_paths = dirname(stan_code_map), 
                                                               cpp_options = list(stan_threads = TRUE, STAN_CPP_OPTIMS = TRUE))
    idealdata@stanmodel_gpu <- stan_code_gpu %>% cmdstan_model(include_paths = dirname(stan_code_map), 
                                                               cpp_options = list(stan_threads = TRUE, STAN_CPP_OPTIMS = TRUE, 
                                                                                  STAN_OPENCL = TRUE, opencl_platform_id = 0, opencl_device_id = 0))
  } else {
    idealdata@stanmodel_map <- stan_code_map %>% cmdstan_model(include_paths = dirname(stan_code_map), 
                                                               cpp_options = list(stan_threads = TRUE))
    idealdata@stanmodel_gpu <- stan_code_gpu %>% cmdstan_model(include_paths = dirname(stan_code_map), 
                                                               cpp_options = list(stan_threads = TRUE, STAN_OPENCL = TRUE, 
                                                                                  opencl_platform_id = 0, opencl_device_id = 0))
  }
  if (vary_ideal_pts == "none") {
    idealdata@score_matrix$time_id <- 1
  }
  vary_ideal_pts <- switch(vary_ideal_pts, none = 1, random_walk = 2, 
                           AR1 = 3, GP = 4)
  mod_count <- length(unique(idealdata@score_matrix$model_id))
  if (length(unique(idealdata@score_matrix$ordered_id)) > 1) {
    mod_count <- mod_count + length(unique(idealdata@score_matrix$ordered_id)) - 
      1
  }
  if (idealdata@score_matrix$model_id[1] == "missing") {
    idealdata@score_matrix$model_id <- model_type
    idealdata@score_matrix$discrete <- as.numeric(model_type %in% 
                                                    c(1, 2, 3, 4, 5, 6, 7, 8, 13, 14))
  } else {
    if (!is.numeric(idealdata@score_matrix$model_id)) {
      stop("Column model_id in the data is not a numeric series of \n           integers. Please pass a numeric value in the id_make function\n           for model_id based on the available model types.")
    }
  }
  if ((!(all(c(restrict_ind_high, restrict_ind_low) %in% unique(idealdata@score_matrix$person_id))) && 
       !use_groups) && const_type == "persons") {
    stop("Your restricted persons/items are not in the data.")
  } else if (!(all(c(restrict_ind_high, restrict_ind_low) %in% 
                 unique(idealdata@score_matrix$item_id))) && const_type == 
           "items") {
    stop("Your restricted persons/items are not in the data.")
  } else if ((!(all(c(restrict_ind_high, restrict_ind_low) %in% 
                  unique(idealdata@score_matrix$group_id))) && use_groups) && 
           const_type == "persons") {
    stop("Your restricted persons/items are not in the data.")
  }
  if (any(c(3, 4, 5, 6) %in% idealdata@score_matrix$model_id)) {
    if (is.null(idealdata@score_matrix$ordered_id) || !is.numeric(idealdata@score_matrix$ordered_id)) {
      stop("If you have an ordered categorical variable in a multi-distribution model, you must include a column in the data with the number of ordered categories for each row in the data.\n\n             See id_make documentation for more info.")
    }
  } else {
    idealdata@score_matrix$ordered_id <- 0
  }
  if (use_groups == T) {
    legispoints <- as.numeric(idealdata@score_matrix$group_id)
  } else {
    legispoints <- as.numeric(idealdata@score_matrix$person_id)
  }
  billpoints <- as.numeric(idealdata@score_matrix$item_id)
  timepoints <- as.numeric(factor(as.numeric(idealdata@score_matrix$time_id)))
  modelpoints <- as.integer(idealdata@score_matrix$model_id)
  ordered_id <- as.integer(idealdata@score_matrix$ordered_id)
  time_ind <- switch(class(idealdata@score_matrix$time_id)[1], 
                     factor = unique(as.numeric(idealdata@score_matrix$time_id)), 
                     Date = unique(as.numeric(idealdata@score_matrix$time_id)), 
                     POSIXct = unique(as.numeric(idealdata@score_matrix$time_id)), 
                     POSIXlt = unique(as.numeric(idealdata@score_matrix$time_id)), 
                     numeric = unique(idealdata@score_matrix$time_id), integer = unique(idealdata@score_matrix$time_id))
  if (gp_min_length >= gp_num_diff[1]) {
    stop("The parameter gp_min_length cannot be equal to or greater than gp_num_diff[1].")
  }
  if (("outcome_cont" %in% names(idealdata@score_matrix)) && 
      ("outcome_disc" %in% names(idealdata@score_matrix))) {
    Y_int <- idealdata@score_matrix$outcome_disc
    Y_cont <- idealdata@score_matrix$outcome_cont
  } else if ("outcome_cont" %in% names(idealdata@score_matrix)) {
    Y_int <- array(0)
    Y_cont <- idealdata@score_matrix$outcome_cont
  } else {
    Y_cont <- array(0)
    Y_int <- idealdata@score_matrix$outcome_disc
  }
  remove_list <- idealstan:::.remove_nas(Y_int, Y_cont, discrete = idealdata@score_matrix$discrete, 
                             legispoints, billpoints, timepoints, modelpoints, ordered_id, 
                             idealdata, time_ind = as.array(time_ind), time_proc = vary_ideal_pts, 
                             ar_sd = ar_sd, gp_sd_par = gp_sd_par, num_diff = gp_num_diff, 
                             m_sd_par = gp_m_sd_par, min_length = gp_min_length, const_type = switch(const_type, 
                                                                                                     persons = 1L, items = 2L), discrim_reg_sd = discrim_reg_sd, 
                             discrim_miss_sd = discrim_miss_sd, diff_reg_sd = diff_reg_sd, 
                             diff_miss_sd = diff_miss_sd, legis_sd = person_sd, restrict_sd_high = restrict_sd_high, 
                             restrict_sd_low = restrict_sd_low, restrict_high = idealdata@restrict_ind_high, 
                             restrict_low = idealdata@restrict_ind_low, fix_high = idealdata@restrict_num_high, 
                             fix_low = idealdata@restrict_num_low)
  out_list <- idealstan:::.make_sum_vals(idealdata@score_matrix, map_over_id, 
                             use_groups = use_groups, remove_nas = remove_list$remove_nas)
  sum_vals <- out_list$sum_vals
  S <- nrow(sum_vals)
  if (het_var) {
    num_var <- length(unique(remove_list$billpoints[remove_list$modelpoints %in% 
                                                      c(9, 10, 11, 12)]))
    mod_items <- tibble(model_id = remove_list$modelpoints, 
                        item_id = remove_list$billpoints) %>% distinct
    mod_items <- mutate(mod_items, cont = model_id %in% c(9, 
                                                          10, 11, 12)) %>% group_by(cont) %>% mutate(num_var = 1:n())
    type_het_var <- arrange(mod_items, item_id) %>% pull(num_var)
  } else {
    num_var <- 1
    type_het_var <- rep(num_var, length(unique(billpoints)))
  }
  if (!is.null(boundary_prior)) {
    if (is.null(boundary_prior$beta)) {
      stop("If you want to use a boundary-avoiding prior for time-series variance, please pass a list with an element named beta, i.e. list(beta=1).")
    }
    if (boundary_prior$beta <= 0) {
      stop("Boundary prior beta value must be strictly positive (i.e. greater than 0).")
    }
    inv_gamma_beta <- boundary_prior$beta
  } else {
    inv_gamma_beta <- 0
  }
  if (remove_list$N_cont > 0) {
    Y_cont <- remove_list$Y_cont[out_list$this_data$orig_order]
  } else {
    Y_cont <- remove_list$Y_cont
  }
  if (remove_list$N_int > 0) {
    Y_int <- remove_list$Y_int[out_list$this_data$orig_order]
  } else {
    Y_int <- remove_list$Y_int
  }
  if (!is.null(ignore_db)) {
    if (!all(c("person_id", "time_id", "ignore") %in% names(ignore_db))) {
      stop("Ignore DB does not have the three necessary columns: time_id, person_id, and ignore (0 or 1).")
    }
    ignore_db <- mutate(ungroup(ignore_db), person_id = as.numeric(person_id), 
                        time_id = as.numeric(factor(as.numeric(ignore_db$time_id))))
    ignore_db <- filter(ignore_db, time_id %in% remove_list$timepoints, 
                        person_id %in% remove_list$legispoints) %>% select(person_id, 
                                                                           time_id, ignore) %>% arrange(person_id, time_id) %>% 
      as.matrix
  } else {
    ignore_db <- matrix(nrow = 0, ncol = 3)
  }
  this_data <- list(N = remove_list$N, N_cont = remove_list$N_cont, 
                    N_int = remove_list$N_int, Y_int = Y_int, Y_cont = Y_cont, 
                    y_int_miss = remove_list$y_int_miss, y_cont_miss = remove_list$y_cont_miss, 
                    num_var = num_var, type_het_var = type_het_var, S = nrow(sum_vals), 
                    S_type = as.numeric(map_over_id == "persons"), T = remove_list$max_t, 
                    num_legis = remove_list$num_legis, num_bills = remove_list$num_bills, 
                    num_ls = remove_list$num_ls, num_bills_grm = remove_list$num_bills_grm, 
                    ll = remove_list$legispoints[out_list$this_data$orig_order], 
                    bb = remove_list$billpoints[out_list$this_data$orig_order], 
                    mm = remove_list$modelpoints[out_list$this_data$orig_order], 
                    ignore = as.numeric(nrow(ignore_db) > 0), ignore_db = ignore_db, 
                    mod_count = length(unique(remove_list$modelpoints)), 
                    num_fix_high = as.integer(1), num_fix_low = as.integer(1), 
                    tot_cats = length(remove_list$n_cats_rat), n_cats_rat = remove_list$n_cats_rat, 
                    n_cats_grm = remove_list$n_cats_grm, order_cats_rat = remove_list$order_cats_rat[out_list$this_data$orig_order], 
                    order_cats_grm = remove_list$order_cats_grm[out_list$this_data$orig_order], 
                    # num_bills_grm = ifelse(any(remove_list$modelpoints %in% 
                    #                              c(5, 6)), remove_list$num_bills, 0L), 
                    LX = remove_list$LX, 
                    SRX = remove_list$SRX, SAX = remove_list$SAX, legis_pred = remove_list$legis_pred[out_list$this_data$orig_order, 
                                                                                                      , drop = FALSE], srx_pred = remove_list$srx_pred[out_list$this_data$orig_order, 
                                                                                                                                                       , drop = FALSE], sax_pred = remove_list$sax_pred[out_list$this_data$orig_order, 
                                                                                                                                                                                                        , drop = FALSE], time = remove_list$timepoints[out_list$this_data$orig_order], 
                    time_proc = vary_ideal_pts, 
                    discrim_reg_sd = discrim_reg_sd, 
                    discrim_abs_sd = discrim_miss_sd, diff_reg_sd = diff_reg_sd, 
                    diff_abs_sd = diff_miss_sd, legis_sd = person_sd, restrict_sd_high = 5, 
                    restrict_sd_low = 5, time_sd = time_fix_sd, time_var_sd = time_var, 
                    ar1_up = ar1_up, ar1_down = ar1_down, inv_gamma_beta = inv_gamma_beta, 
                    center_cutoff = as.integer(time_center_cutoff), restrict_var = restrict_var, 
                    ar_sd = ar_sd, zeroes = as.numeric(inflate_zero), time_ind = as.array(time_ind), 
                    time_proc = vary_ideal_pts, 
                    gp_sd_par = gp_sd_par, m_sd_par = gp_m_sd_par, 
                    min_length = gp_min_length, id_refresh = id_refresh, 
                    sum_vals = as.matrix(sum_vals), const_type = switch(const_type, 
                                                                        persons = 1L, items = 2L), restrict_high = 1, restrict_low = 2, 
                    fix_high = 0, fix_low = 0, num_diff = gp_num_diff, pos_discrim = as.integer(pos_discrim), 
                    grainsize = grainsize)
  idealdata <- idealstan:::id_model(object = idealdata, fixtype = fixtype, 
                        this_data = this_data, nfix = nfix, restrict_ind_high = restrict_ind_high, 
                        restrict_ind_low = restrict_ind_low, ncores = ncores, 
                        tol_rel_obj = tol_rel_obj, use_groups = use_groups, prior_fit = prior_fit, 
                        fix_high = fix_high, fix_low = fix_low, const_type = const_type)
  if (("outcome_cont" %in% names(idealdata@score_matrix)) && 
      ("outcome_disc" %in% names(idealdata@score_matrix))) {
    Y_int <- idealdata@score_matrix$outcome_disc
    Y_cont <- idealdata@score_matrix$outcome_cont
  } else if ("outcome_cont" %in% names(idealdata@score_matrix)) {
    Y_int <- array(0)
    Y_cont <- idealdata@score_matrix$outcome_cont
  } else {
    Y_cont <- array(0)
    Y_int <- idealdata@score_matrix$outcome_disc
  }
  remove_list <- idealstan:::.remove_nas(Y_int, Y_cont, discrete = idealdata@score_matrix$discrete, 
                             legispoints, billpoints, timepoints, modelpoints, ordered_id, 
                             idealdata, time_ind = as.array(time_ind), time_proc = vary_ideal_pts, 
                             ar_sd = ar_sd, gp_sd_par = gp_sd_par, m_sd_par = gp_m_sd_par, 
                             min_length = gp_min_length, const_type = switch(const_type, 
                                                                             persons = 1L, items = 2L), discrim_reg_sd = discrim_reg_sd, 
                             discrim_miss_sd = discrim_miss_sd, diff_reg_sd = diff_reg_sd, 
                             diff_miss_sd = diff_miss_sd, legis_sd = person_sd, restrict_sd_high = restrict_sd_high, 
                             restrict_sd_low = restrict_sd_low, restrict_high = idealdata@restrict_ind_high, 
                             restrict_low = idealdata@restrict_ind_low, fix_high = idealdata@restrict_num_high, 
                             fix_low = idealdata@restrict_num_low)
  idealdata <- remove_list$idealdata
  out_list <- idealstan:::.make_sum_vals(idealdata@score_matrix, map_over_id, 
                             use_groups = use_groups, remove_nas = remove_list$remove_nas)
  sum_vals <- out_list$sum_vals
  idealdata@score_matrix <- out_list$this_data
  S <- nrow(sum_vals)
  if (het_var) {
    num_var <- length(unique(remove_list$billpoints[remove_list$modelpoints %in% 
                                                      c(9, 10, 11, 12)]))
    mod_items <- tibble(model_id = this_data$mm, item_id = this_data$bb) %>% 
      distinct
    mod_items <- mutate(mod_items, cont = model_id %in% c(9, 
                                                          10, 11, 12)) %>% group_by(cont) %>% mutate(num_var = 1:n())
    type_het_var <- arrange(mod_items, item_id) %>% pull(num_var)
  } else {
    num_var <- 1
    type_het_var <- rep(num_var, length(unique(billpoints)))
  }
  if (remove_list$N_cont > 0) {
    Y_cont <- remove_list$Y_cont[out_list$this_data$orig_order]
  } else {
    Y_cont <- remove_list$Y_cont
  }
  if (remove_list$N_int > 0) {
    Y_int <- remove_list$Y_int[out_list$this_data$orig_order]
  } else {
    Y_int <- remove_list$Y_int
  }
  this_data <- list(N = remove_list$N, N_cont = remove_list$N_cont, 
                    N_int = remove_list$N_int, Y_int = Y_int, Y_cont = Y_cont, 
                    y_int_miss = remove_list$y_int_miss, y_cont_miss = remove_list$y_cont_miss, 
                    num_var = num_var, type_het_var = type_het_var, S = nrow(sum_vals), 
                    S_type = as.numeric(map_over_id == "persons"), T = remove_list$max_t, 
                    num_legis = remove_list$num_legis, num_bills = remove_list$num_bills, 
                    num_ls = remove_list$num_ls, 
                    num_bills_grm = remove_list$num_bills_grm, 
                    ll = remove_list$legispoints[out_list$this_data$orig_order], 
                    bb = remove_list$billpoints[out_list$this_data$orig_order], 
                    mm = remove_list$modelpoints[out_list$this_data$orig_order], 
                    ignore = as.numeric(nrow(ignore_db) > 0), ignore_db = ignore_db, 
                    mod_count = length(unique(remove_list$modelpoints)), 
                    num_fix_high = as.integer(1), num_fix_low = as.integer(1), 
                    tot_cats = length(remove_list$n_cats_rat), n_cats_rat = remove_list$n_cats_rat, 
                    n_cats_grm = remove_list$n_cats_grm, order_cats_rat = remove_list$order_cats_rat[out_list$this_data$orig_order], 
                    order_cats_grm = remove_list$order_cats_grm[out_list$this_data$orig_order], 
                    # num_bills_grm = ifelse(any(remove_list$modelpoints %in% 
                    #                              c(5, 6)), remove_list$num_bills, 0L),
                    LX = remove_list$LX, 
                    SRX = remove_list$SRX, SAX = remove_list$SAX, legis_pred = remove_list$legis_pred[out_list$this_data$orig_order, 
                                                                                                      , drop = FALSE], srx_pred = remove_list$srx_pred[out_list$this_data$orig_order, 
                                                                                                                                                       , drop = FALSE], sax_pred = remove_list$sax_pred[out_list$this_data$orig_order, 
                                                                                                                                                                                                        , drop = FALSE], time = remove_list$timepoints[out_list$this_data$orig_order], 
                    # time_proc = vary_ideal_pts, 
                    discrim_reg_sd = discrim_reg_sd, 
                    discrim_abs_sd = discrim_miss_sd, diff_reg_sd = diff_reg_sd, 
                    diff_abs_sd = diff_miss_sd, legis_sd = person_sd, restrict_sd_high = restrict_sd_high, 
                    restrict_sd_low = restrict_sd_low, time_sd = time_fix_sd, 
                    time_var_sd = time_var, ar1_up = ar1_up, ar1_down = ar1_down, 
                    inv_gamma_beta = inv_gamma_beta, center_cutoff = as.integer(time_center_cutoff), 
                    restrict_var = restrict_var, ar_sd = ar_sd, zeroes = as.numeric(inflate_zero), 
                    time_ind = as.array(time_ind), 
                    time_proc = vary_ideal_pts, 
                    gp_sd_par = gp_sd_par, m_sd_par = gp_m_sd_par, min_length = gp_min_length, 
                    id_refresh = id_refresh, sum_vals = as.matrix(sum_vals), 
                    const_type = switch(const_type, persons = 1L, items = 2L), 
                    restrict_high = idealdata@restrict_ind_high, restrict_low = idealdata@restrict_ind_low, 
                    fix_high = idealdata@restrict_num_high, fix_low = idealdata@restrict_num_low, 
                    num_diff = gp_num_diff, pos_discrim = as.integer(pos_discrim), 
                    grainsize = grainsize)
  idealdata@n_cats_rat <- remove_list$n_cats_rat
  idealdata@n_cats_grm <- remove_list$n_cats_grm
  idealdata@order_cats_rat <- remove_list$order_cats_rat
  idealdata@order_cats_grm <- remove_list$order_cats_grm
  outobj <- idealstan:::sample_model(object = idealdata, nchains = nchains, 
                         niters = niters, warmup = warmup, ncores = ncores, this_data = this_data, 
                         use_vb = use_vb, gpu = gpu, save_files = save_files, 
                         keep_param = keep_param, tol_rel_obj = tol_rel_obj, within_chain = within_chain)
  outobj@model_type <- model_type
  outobj@time_proc <- vary_ideal_pts
  outobj@use_groups <- use_groups
  outobj@map_over_id <- map_over_id
  outobj@time_fix_sd <- time_fix_sd
  outobj@restrict_var <- restrict_var
  outobj@time_center_cutoff <- time_center_cutoff
  if (this_data$T > 1 && ((!is.null(keep_param$person_vary) && 
                           keep_param$person_vary) || is.null(keep_param))) {
    outobj@time_varying <- try(idealstan:::.get_varying(outobj, legis_x = remove_list$legis_pred[out_list$this_data$orig_order, 
                                                                                     , drop = FALSE], person_id = this_data$ll, time_id = this_data$time))
  }
  return(outobj)
}

# Function to extract ideal point values

extract_ideal_points <- function (object, return_data = FALSE, include = NULL, item_plot = NULL, 
                                  text_size_label = 2, text_size_group = 2.5, high_limit = 0.95, 
                                  low_limit = 0.05, line_size = 1, highlight = NULL, plot_text = TRUE, 
                                  use_ci = TRUE, person_line_alpha = 0.3, person_ci_alpha = 0.8, 
                                  item_plot_type = "non-inflated", show_true = FALSE, group_color = TRUE, 
                                  hpd_limit = 10, sample_persons = NULL, plot_sim = FALSE, 
                                  ...) 
{
  person_labels <- quo(person_id)
  group_labels <- quo(group_id)
  already_facet <- FALSE
  if (class(object) == "idealstan") {
    person_params <- prepare_legis_data(object, high_limit = high_limit, 
                                        low_limit = low_limit)
    model_wrap <- FALSE
    use_groups <- object@use_groups
  }
  else {
    person_params <- lapply(object, .prepare_legis_data, 
                            high_limit = high_limit, low_limit = low_limit) %>% 
      bind_rows(.id = "Model")
    check_groups <- sapply(object, function(x) x@use_groups)
    if (all(check_groups)) {
      use_groups <- T
    }
    else {
      use_groups <- F
    }
    model_wrap <- TRUE
    object <- object[[1]]
  }
  if (!is.null(item_plot)) {
    if (object@model_type %in% c(1, 2) || (object@model_type > 
                                           6 && object@model_type < 13)) {
      item_points <- lapply(item_plot, .item_plot_binary, 
                            object = object, low_limit = low_limit, high_limit = high_limit) %>% 
        bind_rows()
    }
    else if (object@model_type %in% c(3, 4)) {
      item_points <- lapply(item_plot, .item_plot_ord_rs, 
                            object = object, low_limit = low_limit, high_limit = high_limit) %>% 
        bind_rows()
    }
    else if (object@model_type %in% c(5, 6)) {
      item_points <- lapply(item_plot, .item_plot_ord_grm, 
                            object = object, low_limit = low_limit, high_limit = high_limit) %>% 
        bind_rows()
    }
    else if (object@model_type %in% c(13, 14)) {
      item_points <- lapply(item_plot, .item_plot_ls, object = object, 
                            low_limit = low_limit, high_limit = high_limit) %>% 
        bind_rows()
    }
    item_points <- left_join(item_points, object@score_data@score_matrix, 
                             by = c(item_name = "item_id")) %>% gather(key = "ci_type", 
                                                                       value = "ci_value", item_high, item_low) %>% mutate(ci_type = factor(ci_type, 
                                                                                                                                            levels = c("item_low", "item_high"), labels = c(paste0(low_limit * 
                                                                                                                                                                                                     100, "%"), paste0(high_limit * 100, "%"))))
    person_params <- left_join(person_params, item_points, 
                               by = c(person_id = "person_id", group_id = "group_id"))
    person_params$time_id <- person_params$time_id.x
    person_params$time_id.x <- NULL
    person_params$time_id.y <- NULL
    person_params <- spread(person_params, key = ci_type, 
                            value = ci_value)
    if (item_plot_type != "both") {
      item_plot_type <- switch(item_plot_type, `non-inflated` = "Non-Inflated\nDiscrimination", 
                               inflated = "Inflated\nDiscrimination")
      person_params <- filter(person_params, item_type == 
                                item_plot_type)
    }
  }
  if (!is.null(include)) {
    if (use_groups) {
      person_params <- filter(person_params, group_id %in% 
                                include)
    }
    else {
      person_params <- filter(person_params, person_id %in% 
                                include)
    }
  }
  if (use_groups && !model_wrap) {
    base_id <- ~group_id
  }
  else if (!use_groups && !model_wrap) {
    base_id <- ~person_id
  }
  else {
    base_id <- ~Model
    if (use_groups) {
      wrap_id <- ~group_id
    }
    else {
      wrap_id <- ~person_id
    }
  }
  if (!is.null(object@score_data@simul_data) && plot_sim == 
      T) {
    true_pts <- object@score_data@simul_data$true_person
    colnames(true_pts) <- c(as.character(1:ncol(true_pts)))
    true_pts <- as_data_frame(true_pts) %>% mutate(person_id = 1:n()) %>% 
      gather(key = time_id, value = true_pt, -person_id) %>% 
      mutate(time_id = as.numeric(time_id), person_id = factor(person_id), 
             person_id = fct_relevel(person_id, as.character(object@score_data@restrict_ind_low), 
                                     as.character(object@score_data@restrict_ind_high), 
                                     after = length(levels(person_id))))
    person_params <- left_join(person_params, true_pts, by = c("person_id", 
                                                               "time_id"))
  }
  if (any(person_params$high_pt > 1e+06)) {
    warnings("Very large person ideal points detected. Removing from plot, likely a sign of model misfit.")
    person_params <- filter(person_params, high_pt < 1e+06, 
                            low_pt > (-1e+06))
  }
  if (use_ci == T) {
    outplot <- person_params %>% ggplot(aes_(x = ~time_id)) + 
      geom_ribbon(aes_(ymin = ~low_pt, ymax = ~high_pt, 
                       group = base_id), fill = "grey80", colour = NA, 
                  alpha = person_ci_alpha)
  }
  else {
    outplot <- person_params %>% ggplot(aes_(x = ~time_id))
  }
  if (!is.null(object@score_data@simul_data) && plot_sim == 
      T) {
    outplot <- outplot + geom_line(aes_(y = ~true_pt, colour = base_id), 
                                   alpha = person_ci_alpha, size = line_size)
  }
  else {
    if (group_color) {
      if (model_wrap) {
        outplot <- outplot + geom_line(aes_(y = ~median_pt, 
                                            group = base_id, colour = ~Model), alpha = person_ci_alpha, 
                                       size = line_size)
      }
      else {
        outplot <- outplot + geom_line(aes_(y = ~median_pt, 
                                            group = base_id, colour = ~group_id), alpha = person_ci_alpha, 
                                       size = line_size)
      }
    }
    else {
      outplot <- outplot + geom_line(aes_(y = ~median_pt, 
                                          group = base_id), alpha = person_ci_alpha, size = line_size)
    }
  }
  if (!is.null(item_plot)) {
    outplot <- outplot + geom_hline(aes(yintercept = item_median), 
                                    linetype = 2, alpha = 0.6) + geom_ribbon(aes(ymin = `5%`, 
                                                                                 ymax = `95%`), alpha = 0.2, fill = "pink") + guides(colour = guide_legend("")) + 
      theme(legend.position = "bottom")
    if (item_plot_type == "both" && length(item_plot) > 1) {
      if (model_wrap) {
        outplot <- outplot + facet_wrap(update(wrap_id, 
                                               ~. + item_name + item_type), dir = "v", ncol = length(pull(person_params, 
                                                                                                          Model)))
      }
      else {
        outplot <- outplot + facet_wrap(~item_name + 
                                          item_type, dir = "v")
      }
    }
    else if (item_plot_type == "both") {
      if (model_wrap) {
        outplot <- outplot + facet_wrap(update(wrap_id, 
                                               ~. + item_type), dir = "v", ncol = length(pull(person_params, 
                                                                                              Model)))
      }
      else {
        outplot <- outplot + facet_wrap(~item_type, dir = "v")
      }
    }
    already_facet <- TRUE
  }
  if (plot_text == TRUE) {
    if (model_wrap) {
      sampled_data <- group_by(person_params, !!as_quosure(base_id), 
                               !!as_quosure(wrap_id)) %>% sample_n(1)
    }
    else {
      sampled_data <- group_by(person_params, !!as_quosure(base_id)) %>% 
        sample_n(1)
    }
    if (!is.null(highlight)) {
      sampled_data <- filter(sampled_data, !(!!as_quosure(base_id) %in% 
                                               highlight))
    }
    outplot <- outplot + geom_text_repel(aes_(x = ~time_id, 
                                              y = ~median_pt, label = base_id), data = sampled_data, 
                                         size = text_size_label, segment.colour = "grey50", 
                                         segment.alpha = 0.5)
  }
  if (!is.null(highlight)) {
    outplot <- outplot + gghighlight(!!as_quosure(base_id) %in% 
                                       highlight, use_group_by = F)
  }
  if (group_color == F || (group_color == F && is.null(highlight))) {
    output <- outplot + guides(color = FALSE, fill = FALSE)
  }
  outplot <- outplot + theme_minimal() + ylab("Ideal Point Scale") + 
    xlab("") + theme(panel.grid = element_blank())
  if (!already_facet && model_wrap) {
    outplot <- outplot + facet_wrap(wrap_id, scales = "free_y", 
                                    ncol = length(pull(person_params, Model)))
  }
  return(outplot)
}

# Function to tidy data from ideal points estimation

prepare_legis_data <- 
  function (object, high_limit = NULL, low_limit = NULL, aggregate = TRUE, 
            type = "ideal_pts") 
  {
    if (length(unique(object@score_data@score_matrix$time_id)) > 
        1 && type != "variance") {
      person_params <- object@time_varying
      if ("draws_matrix" %in% class(person_params)) {
        person_params <- person_params %>% as_draws_df %>% 
          dplyr::select(-.chain, -.iteration, -.draw)
      }
      else {
        person_params <- as_tibble(person_params)
      }
      if (aggregate) {
        person_params <- person_params %>% 
          gather(key = legis, value = ideal_pts) %>% 
          group_by(legis) %>% 
          summarize(low_pt = quantile(ideal_pts, low_limit, na.rm = T), 
                    high_pt = quantile(ideal_pts, high_limit, na.rm = T), 
                    median_pt = median(ideal_pts, na.rm = T)) %>% 
          mutate(param_id = stringr::str_extract(legis, "[0-9]+\\]"), 
                 param_id = as.numeric(stringr::str_extract(param_id, "[0-9]+")), 
                 time_point = stringr::str_extract(legis, "\\[[0-9]+"), 
                 time_point = as.numeric(stringr::str_extract(time_point, "[0-9]+")))
      }
      else {
        person_params <- person_params %>% gather(key = legis, 
                                                  value = ideal_pts) %>% group_by(legis) %>% mutate(param_id = stringr::str_extract(legis, 
                                                                                                                                    "[0-9]+\\]"), param_id = as.numeric(stringr::str_extract(param_id, 
                                                                                                                                                                                             "[0-9]+")), time_point = stringr::str_extract(legis, 
                                                                                                                                                                                                                                           "\\[[0-9]+"), time_point = as.numeric(stringr::str_extract(time_point, 
                                                                                                                                                                                                                                                                                                      "[0-9]+")))
      }
      person_ids <- select(object@score_data@score_matrix, 
                           !!quo(person_id), !!quo(time_id), !!quo(group_id)) %>% 
        distinct %>% mutate(person_id_num = as.numeric(!!quo(person_id)), 
                            time_id_num = as.numeric(factor(!!quo(time_id))), 
                            group_id_num = as.numeric(!!quo(group_id)))
      if (object@use_groups) {
        person_params <- person_params %>% left_join(person_ids, 
                                                     by = c(param_id = "group_id_num", time_point = "time_id_num"))
        person_params <- person_params %>% group_by(param_id) %>% 
          arrange(param_id, time_point) %>% mutate(time_id = rep(sort(unique(person_ids$time_id)), 
                                                                 each = length(unique(person_ids$person_id[person_ids$group_id_num == 
                                                                                                             param_id])))) %>% fill(everything(), .direction = "downup") %>% 
          ungroup
      }
      else {
        person_params <- person_params %>% left_join(person_ids, 
                                                     by = c(param_id = "person_id_num", time_point = "time_id_num"))
        person_params <- person_params %>% group_by(param_id) %>% 
          arrange(param_id, time_point) %>% mutate(time_id = sort(unique(person_ids$time_id))) %>% 
          fill(everything(), .direction = "downup") %>% 
          ungroup
      }
    }
    else {
      if (type == "ideal_pts") {
        person_params <- object@stan_samples$draws("L_full") %>% 
          as_draws_df %>% dplyr::select(-.chain, -.iteration, 
                                        -.draw)
      }
      else if (type == "variance") {
        person_params <- object@stan_samples$draws("time_var_free") %>% 
          as_draws_df %>% dplyr::select(-.chain, -.iteration, 
                                        -.draw)
      }
      person_params <- person_params %>% gather(key = legis, 
                                                value = ideal_pts)
      person_ids <- tibble(long_name = person_params$legis) %>% 
        distinct
      legis_nums <- stringr::str_extract_all(person_ids$long_name, 
                                             "[0-9]+", simplify = T)
      person_ids <- mutate(person_ids, id_num = as.numeric(legis_nums))
      person_data <- distinct(select(object@score_data@score_matrix, 
                                     person_id, group_id))
      if (object@use_groups) {
        person_data <- mutate(person_data, id_num = as.numeric(group_id))
      }
      else {
        person_data <- mutate(person_data, id_num = as.numeric(person_id))
      }
      person_ids <- left_join(person_ids, person_data)
      if (aggregate) {
        person_params <- person_params %>% group_by(legis) %>% 
          summarize(low_pt = quantile(ideal_pts, low_limit), 
                    high_pt = quantile(ideal_pts, high_limit), 
                    median_pt = median(ideal_pts)) %>% left_join(person_ids, 
                                                                 by = c(legis = "long_name"))
      }
      else {
        person_params <- person_params %>% left_join(person_ids, 
                                                     by = c(legis = "long_name"))
      }
    }
    person_params
  }

# Function to median centering

rescale_center_median <- function(my_numbers) {
  my_median = median(my_numbers, na.rm = TRUE)
  my_range = range(my_numbers, na.rm = TRUE)
  scale_factor = max(abs(my_range-my_median))
  (my_numbers - my_median) / scale_factor
}